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1. Introduction 

Understanding the statistical properties of open, many-particles system is one of the 
challenges of nonequilibrium statistical mechanics. From a fundamental point of view, a 
successful approach would require to find, and possibly compute explicitely, a statistical 
measure for (at least) systems steadily kept out of equilibrium. Some insight has been 
gained over the years mostly thanks to the analysis of specific models (for a recent 
account see e.g. [1] and references therein). A related open problem is the derivation 
of phenomenological transport laws from the microscopic dynamics, without any ad hoc 
statistical assumption. An example is the famous law, postulated by Joseph Fourier 
almost two hundred years ago, relating the heat flux J flowing within a solid material 
to the local temperature gradient, 

J=-kVT, (1) 

where the constant of proportionality k, is the thermal conductivity. 
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In the lack of a general framework, simple models are precious to attack such 
difficult problems j2l [3]. An instance, dating back to 1967, was provided by Rieder, 
Lebowitz and Lieb who considered heat conduction in a chain of harmonic oscillators 
connected at its boundaries to two stochastic heat baths [I]. They showed that the 
invariant measure in phase-space (i.e. the stationary solution of the associated Fokker- 
Planck equation) is a multivariate Gaussian. Furthermore, they proved that, due to 
the integrability of the underlying dynamics, such a model is not able to support a 
temperature gradient. However, this is one of the very few systems that have been 
rigorously solved. Extensions of this model, where anharmonicities are introduced 
by means of self-consistent local thermostats, were extensively studied [SI El [7j. In 
recent years, further attempts to derive Fourier's law in deterministic systems have 
been reported [HI El EH H] . 

As a complementary approach, stochastic models have played an important role 
in understanding how energy is microscopically transported. This is mainly due to 
the fact that the stochastic approach seems to easily yield results that would require 
much more efforts by adopting the dynamical approach. In fact, while stochastic 
models are assumed to be a reduced (mesoscopic) representation of the "chaotic" 
microscopic dynamics, they are free from the intricacies of the fractal structures 
arising in deterministic dynamics. Actually, the leap from such class of models to 
even the simplest deterministic, nonlinear ones is still a challenge for the theory 
[3]. At the simplest level of modeling, energy is assumed to be randomly exchanged 
between neighbouring sites of a lattice [12, [131 El US]- This class of systems has 
the invaluable advantage of allowing for a mathematically rigorous treatment, which 
is usually unfeasible in the deterministic case. Recently, systems of harmonic oscillators 
exchanging energy with "conservative" noise have been proven to admit a unique 
stationary state consistent with ([1]) [16] . However, if the additional constraint that 
the random process conserves also linear momentum is imposed, the equilibrium energy- 
current correlation function decays as t~ d l 2 (d being the lattice dimension) and transport 
becomes anomalous in d < 2 [17] . This means that (CD) breaks down as k diverges with 
the system size. The results of [17] thus provide a rigorous basis to the numerical 
evidence of anomalous transport and diffusion in deterministic nonlinear models with 
momentum conservation [2], with the only exception of the coupled rotor chain [T8l fT9]. 

In this paper we consider the problem of heat conduction in a chain of harmonic 
oscillators, coupled at its boundaries with two stochastic heat baths. In addition to 
the deterministic bulk dynamics, we consider a stochastic "noise" dynamical term, 
consisting of collisions occurring at a given rate 7, that exchange the momenta of a 
stochastically chosen pair of neighbour oscillators. The stochastic contribution to the 
dynamics maintains the linearity of the associated Fokker-Planck equation. 

Recently, following a principal component analysis, we have numerically found that, 
in the basis identified by the eigenvectors of the covariance matrix, the nonequilibrium 
invariant measure of this model can be effectively expressed as the product of 
independent distributions aligned along collective modes that are spatially localized with 



Stochastic model of anomalous heat transport 



3 



power-law tails [2U]. Moreover, several variables, such as the amplitudes of these modes, 
turn out to be Gaussian distributed. Accordingly, it appears that the unavoidable 
deviations from a Gaussian behaviour are confined to not-so-relevant observables, 
so that the nonequilibrium invariant measure can be effectively considered to be a 
multivariate Gaussian. Within this aproximation, the covariance matrix provides a 
complete description of the corresponding invariant measure. Here, with the help of a 
suitable continuum limit, we derive leading order expressions for the covariance matrix 
in the steady nonequilibrium state, from which explicit formulae for the temperature 
profile and the energy current, are obtained. It should be noticed that this is the first 
example of an analytic expression for the temperature profile in a system characterized 
by anomalous heat transport (i.e. diverging conductivity). 

The paper is organized as follows. In Section [21 we introduce the stochastic model. 
In Section [3], we define the covariance matrix C and write down the coupled equations 
governing the evolution of C towards its stationary value. The key results of the paper 
are also summarized there. Section H] contains the details of the analytical calculation of 
the stationary covariance matrix in the thermodynamic limit N — > oo. Our approach is 
based on a suitable continuum approximation, in which the finite-difference equations 
for the entries of C are replaced by partial differential equations for the corresponding 
field-like variables. We obtain the covariance matrix to leading order in the smallness 
parameter e — 1/ y/~N. In Section [5] we discuss the physical meaning of the analytic 
expressions, compare them with the numerical solution for finite size chains and briefly 
comment on the open problems. 



2. Stochastic Model 



We consider a homogeneous harmonic chain of ./V oscillators of unit mass and frequency 
lj, in contact with two different stochastic Langevin heat baths at its extrema and 
fixed boundary conditions. The dynamics in the bulk of the chain is governed by the 
Hamiltonian 

N r 9 9 



H{q,p,t) = 



2 



to 



(2) 



Furthermore, the 1-st and iV-th oscillators are coupled to Langevin heat baths at 
temperatures T± = T± AT/2 respectively (T is the average temperature (T + + T_)/2). 
Then the equations of motion become 



(3) 



Qn = Vn 

fin = u 2 (q n+1 - 2q n + g re _i) + £„,i(f+ - \q\) + 8 n>N (£_ - \q N ) , 

where £_ and £ + are independent Wiener processes with zero mean and variance 
2\ksT_ and 2AfcfiT + respectively. The fixed boundary conditions are enforced by setting 
Qo — Qn+i = 0. In addition, the chain undergoes random binary collisions, at a rate 7, 
in which the momenta of a couple of neighbouring oscillators are exchanged. Thus, the 
resulting dynamics conserves both total momentum and energy. 



Stochastic model of anomalous heat transport 



4 



The phase space probability density P(q,p,t) of this model, is a solution of the 
Fokker-Planck equation 

^ = (C + C coU )P. (4) 

The first term describing the evolution of the system, as defined by ([3]) , can be written 

as 

£ P - V (A dX > P + D * \ (5) 

L P-^ [A lJ —— + — — — , (5) 



i.j 



where the 2N vector x = q 2 , . . . , qN,Pi,P2, ■ ■ ■ ,Pn), and the 2N x 2N matrices A 
and D are 

A= (w 2 G Ar) ! D= (° 2AA; B T(R + r/S) ) (6) 
with and 1 the null and unit N x N matrices respectively, 

Rjj = + Si,N) 5 = Sij (Si t i — 6i t N) , (7) 

and G is the negative of the Laplacian, 

= 25ij — Si + i t j — . (8) 

Moreover, we introduce the normalized bath temperatures difference rj = AT/T = 
(T + — T_)/T. Finally, the second term in (jlj) associated to stochastic collisions reads 

JV-l 

CcoiiP = 1^2 • • 'ft'+i'Pi' ...)-P(... ,Pj,Pj+i t •■■)]■ ( 9 ) 
j'=i 

Each term in the sum expresses the probability balance for each elementary process in 
which momenta of the pair j, j + 1 are exchanged. 

As we mentioned above, this type of dynamics with conservative noise, was 
originally introduced in [T7] where, however, only the equilibrium case was studied. 
Here we consider the nonequilibrium situation. Moreover, the collision term (Q we 
consider here has two main differences: first, in the present case, only collisions of pairs 
(instead of triplets) are necessary. Second, in [16J, each evolution step is an infinitesimal 
variation of the momenta onto the constant-energy hypersurface. This allows to define a 
generator of the process as a differential operator acting on the p-space. On the contrary, 
in the present case the process remains intrinsically discontinuous. 

3. Covariance Matrix 

Consider the covariance matrix C for the dynamics ([3]), which we write as 

c=(Z? r ), (10) 




where, 

Uij = (qiqj) , Vij = (piPj) , = (qiPj) , (11) 
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(15) 



are three NxN matrices, the square brackets (.) denote an average over P(q,p,t), 
and f denotes the transpose operation. There is no need to include mean values, since 
(Pi) = (ft) = o. 

Note that the matrices U and V are symmetric by definition. The evolution 
equation for C can be written 

C = Co + Ccou , (12) 

where (see e.g. equation (63) in [2]), 

C = D - AC - CA T . (13) 

The collision term C co « is evaluated upon multiplying (jHJ) by XiXj and thereby integrating 
over phase space. We obtain 

C co « = - 7 ( G zt Z £) (14) 

where the auxiliary NxN matrix W is defined by 

Vi-ij-i + Vi+ij+i - 2Vij i = j 

Vj-ij + Vij+i - 2V i>3 - i-j = -l 

V,,,., • V,. ; ; 2V, ; / ./ 1 

Vi+ij- + Vi_i j + Vij_i + V i>3 - +1 - 4Vjj |z - j | > 1 

Equation (fT2l) is thus exact and closed and describes the approach to the 
nonequilibrium steady state. In the present work we aim at finding its stationary 
solution, which amounts to solving the set of linear equations 

Z f = -Z , (16a) 
V = cu 2 UG + AZR + 7 ZG , (166) 
u 2 (GZ + Z T G) + A (RV + VR) + 7 W = 2Xk B T (R + 77S) . (16c) 

Note that for T + = T~ = T, namely 77 = 0, these equations admit the equilibrium 
solution 

7 rji 

U cq = ^G" 1 , V cq = k B Tl , Z cq = . (17) 

U! 

For 77 7^ 0, analogously to what found in purely stochastic models [13 |2TJ, we expect 
the onset of a non-zero heat flux to be accompanied by the appearance of non-diagonal 
terms. 

In the next Section we solve analytically the problem (H6aH16c|) by means of a 
suitable continuum approximation. The idea is to replace the finite-difference equations 
f U6aH16cj) with a set of partial differential equations. Before entering the technical 
details, it is useful to briefly anticipate the main outcomes of our calculation. The 
temperature field Tj = (pf) along the chain, as a function of the scaled variable 
y = 2i/N — 1 can be expressed as 

T(y) = T + AT Q(y) , (18) 
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where Q(y) is defined through the following series, 

Q(y) = -= — Y n - 3 / 2 cos(—(y + l)) , (19) 

(V8-l)C(3/2) o ^ n ^ 2 } ) 

where C(3/2) = 2.612375348... is the Riemann ^-function. It can be seen that Q(y) 
is an odd function of y such that 0(±1) = =Fl/2. The leading term of the stationary 
energy current (see below for the exact definition) is 



s/N 8( v / 8-l)C(3/2) \j lN 
As a consequence, the effective conductivity is 



J 1 /vrVW 

K= AT/A ~ 8( v / 8- l)C(3/2) V 7 
Comments on the physical meaning of these formulae will be given in the last Section. 



4. Analytical solution 

The solution of fll6aH16cj) can be efficiently determined numerically by exploiting the 
sparsity of the corresponding linear problem, as well as the symmetries of the unknowns 
U, V and Z. This approach has been followed in [20J. Here, we solve the problem 
analytically treating the "lattice" equations in the continuum approximation. It must 
be first recognized that the correct scaling is not known a priori, but rather inferred 
from the numerical solution. Therefore, the correctness of the results has to be checked 
a posteriori by consistency. 



4-1. The continuum limit 

The first step consists in mapping the discrete variables % and j into two suitable 
continuous variables x and y, so that an N x N matrix My can be transformed into 
a "field variable" M.(x, y) and the associated discrete equation turned into a partial 
differential equation. In order to do so, it is first necessary to introduce a smallness 
parameter that vanishes as N — >• oo. In [20], it was found that while neighbouring 
elements along the diagonal differ by 0(1/ N), across the diagonal the difference is 
0(1/ \/N). This suggests defining the smallness parameter as e = l/y/N. In addition, 
it is convenient to introduce a further stretching of the longitudinal variable y so as to 
ensure a constant elongation in the (x, y) representation. This is achieved through the 
following transformation 

{i+])e 2 - 1 
1 - \i-j\e' 1 

that is schematically represented also in figure [TJ The nonlinear transformation 
complicates the expansions along y, but is essential to set the boundary conditions 
correctly. Although ([22"]) is singular for \i — j\ = N, this is harmless, since its location 



* = v= V i 7ri ■ (22) 
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Figure 1. Schematic representation of the mapping from the matrix indexes i,j (left) 
to the continuous variables (x,y) (right). The latter vary in the domain T>, definition 
(|23|) (shaded region). The square symbols represent matrix elements, (bolded along 
the diagonal, i = j). Diagonals are parametrically obtained as x — constant, and 
antidiagonals as y = constant. The denominator in the definition of y (|22[) . takes into 
account that the length of the diagonals depend on their value of x so that, the domain 
of y is independent of x. 



diverges to infinity and is thus located outside the region of interest. In the infinite 
volume limit, the variables (x,y) belong in the domain 

V = {(x,y)\xe[0,oo); ye [-1,1]} (23) 

Note that x = const corresponds to moving along a diagonal direction, x = 
corresponding to the main diagonal. 

In order to determine the continuum limit of the equations (U6aH16cl) . it is necessary 
deal with the infinitesimal changes of x and y that arise from Ai and Aj shifts of % and 
j. It is convenient to introduce the integer functions /, s : Z 2 i— > Z 

f{Ai, Aj) = Ai - Aj , s{Ai, Aj) =At + Aj . (24) 

With the help of these shift functions, the coordinates of a point shifted by (Ai, Aj) 
reads 

/ r i (i + j)e 2 - 1 + se 2 , or , 

x' = x + fe; y'= ) // 2 — 2 25 

1 - (i - j)e 2 - fe 2 

where we assume that % > j to get rid of the absolute value. Accordingly, 

V ' = ( y + T~^) 1 - - £X ) ' (26) 
and, up to fourth order in e, 

y > =( y + e 2 s(l +ex + e 2 x 2 )) (l + e 2 f(l + ex + e 2 x 2 ) + e A f 2 ) , (27) 
which is conveniently written as 

y' = y + e 2 (l + ex + e 2 (x 2 + /)) (s + fy) = y + e 2 R Ls , (28) 
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where 

R u = [l + ex + e 2 (x 2 + /)] (fy + s) . (29) 

With these definitions, an infinitesimal change in x involves terms of 0(e) and an 
infinitesimal change in y generates terms of 0(e 2 ), 0(e 3 ) and 0(e A ). However, for the 
estimate of the leading contributions, it is sufficient to consider Rf s = (1 + ex) (fy + s). 

Altogether, the above relations provide a useful tool for investigating the continuum 
limit. For later applications, the above results are summarized in the rule, 

M i+Aitj+Aj = M{x + fe,y + e 2 R f , s ). (30) 

that is written in a convenient form for an expansion in powers of e. Here and in what 
follows, we keep the bold-face notation for the continuous functions derived from the 
matrix variables. 



4-2. Field variables 

The disadvantage of representation (fTTj) is that U cq is a full matrix whose diagonal 
elements are O(N). This hinders the formulation of a proper perturbation scheme 
to compute the non-equilibrium corrections. For the sake of the numerical solution 
carried out in [20], this difficulty has been overcome by looking at correlators involving 
relative rather than absolute displacements, i.e., Z£ • = ((<& — q i+1 )pj) and = 
{(qi+i — qi)(qj+i — qj)). In fact, in this representation, U' turns out to be diagonal 
at equilibrium with diagonal elements of 0(1). On the other hand, Z[j loses the 
antisymmetry of Zjj, a very useful property for our analytical treatment. Therefore, 
we have decided to keep the definition of Z as in ( TTTT) while introducing a new matrix 
Yjj = ui 2 {(q i+ i — <7j)(<7j + i — qj)), which is conveniently expressed in terms of U as, 

Yjj = u; 2 [Ujj — Ujj+i — Uj+ij + Uj+ij+i] . (31) 

Note that the diagonal elements are proportional to the average bond potential energy 

Y i;i = u 2 {{q i+1 - qi ) 2 ) =2^ . (32) 

The next step consists in choosing the proper order of magnitude of the three fields 
V, Y, and Z. This will be done by exploiting the knowledge of the equilibrium case 
and the information arising from previous numerical solution [20]. First of all, since 
Yj j and V i;i are proportional to the mean potential and kinetic energy, respectively, 
they are both of 0(1) as in equilibrium. On the other hand, the off-diagonal elements 
turn out to be of 0(e). Hence, for consistency of the continuum approximation, we 
must consider independently diagonal and off-diagonal (bulk) entries of V and Y. The 
matrix Z exhibits a somehow complementary behaviour. Since it is antisymmetric in x, 
there are no diagonal terms, 

Z(0,y) = 0, (33) 
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while the numerics suggests that in the bulk it is 0(1). We thus define the following 
field variables: in the bulk (i ^ j, x > 0) 

Y id = eY{x, y) + h.o.t. , V<j = eV{x, y) + h.o.t. , = Z{x, y) + h.o.t. . (34) 

and for the diagonal {i = j, x = 0) terms 

V M = T(y) + h.o.t. , Y i(i = 2${y) + h.o.t. , (35) 

The scaling properties of the first corrections to the leading order are not known and the 
comparison with the numerical results discussed in the final section show the existence 
of a singular dependence on e. As a consequence, it is not possible to set up a standard 
perturbation expansion scheme and it is therefore necessary to rely on expressions 
dominated by the leading contributions. In the next section we demonstrate that by 
manipulating equations (I16ffllll66lll6cp and suitable combinations of them, it is possible to 
obtain a set of partial differential equations whose solution gives the covariance matrices 
at leading order. 

4-3. Stationary equation in the bulk 

Using the definitions (jSJ), © and (jHJ), the equation (1 1 6 £>[) writes 

u 2 (2Uij - Uij+x - Uij-i) - Vy + 7 (2Z M - - Zij_i - Z i)i+1 ) = . (36) 
The reader should note that the term ZR appearing in (11661) . can be written as 

AZR= A(Z(-z,-l) + Z(x,l)) . (37) 

This term does only contribute at the boundaries and consequently, we have omitted it 
in (1361) . In the subsequent treatment, this omission will be justified later when we fix 
the boundary conditions of Z. 

In order to write the equations in terms of the new variable Y let us first rewrite 
(I16fej) with i replaced by i + 1: 

u 2 (2Uj + ij — Uj+ij+i — Uj + ij_i) — Vi+ij + 7 (2Zj + ij — Zj+ij-i — Zj + ij + i) = . (38) 
Subtracting ( 1381) from ( 136|) . and using the definition of the matrix Y we obtain 
Yjj — Yjj_i + Vj + ij — Vjj + 7 [— 2Zj + ij + 2Zjj + Zj + ij_i — Zjj_i+ 



Zi+ij+i — — . 



(39) 

With the help of rule ( l30l . the continuous version of (1391) in the bulk is readily written 

as 

Y(x, y) -Y(x + e,y + e 2 i?i,_i) + V(x + e, y + e 2 J? 1;1 ) - V(x, y) + 

7 [-2Z(x + e, y + e 2 i? l5 + 2Z(x, y) + Z(x + 2e,y + e 2 R 2fi )- (40) 

Z(x + e, y + e^i.-i) + Z(x, y + £ 2 -R ,2) -Z(x-e,y + £ 2 R~i,i)] = . 

The leading order of ( 1401) is of 0(e 2 ) and can be written as 

n x (x,y) = 0, (41) 
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where the subscripts denote the partial derivative with respect to the corresponding 
variable and for reasons that will be clear below, we have introduced the function 

n(x,y)=Y(x,y)-V(x,y) . (42) 

Furthermore, by exchanging i with j in equation (139!) and adding the result to (139|) . we 
find a symmetrized equation in the bulk, given by 

2Y i,j - Y m-i - Y i-i,i + + - 2Vjj + 7 [Zjj+i - Z mj + 

Zj+ij-i — Zj_!j +1 + Zj_ij — Zjj_i] = . 



The continuous version of (143]) is 

2Y(x, y) - Y(x + 5, y + e 2 R h _ 1 ) -Y(x-e,y + e 2 R- 1 ,- 1 ) - 2V(x, y) + 

V(x + e,y + £ 2 Ri,x) + Y(x - e, y + e 2 i?_i,i) + 7 [Z(x - e, 3/ + £ 2 i?-i,i)- ^ 

Z(x + e, y + e 2 i?i,i) + Z(x + 2e,y + e 2 R 2fi ) - Z(x -2e,y + e 2 R. 2 ,o)+ 

Z(x -e,y + e 2 R- 1 - 1 ) - Z(x + e,y + e 2 Ri _i)] = , 

whose leading contribution yields 

- fl xx (x, y) + 2 [Y y (x, y) + V y (z, y)] + 2 7 Z xxx (x, y) = . (45) 

By using (HTl) . the above equation becomes 

Y y (x, y ) + V y (x, y ) + 7 Z XXX (x, y ) = . (46) 

Furthermore, integrating (14~T1) on x we obtain that f2(x, y) does not depend on the 
transversal coordinate x, namely 

n(x,y) = ^(y). (47) 

By using this in (1161) to replace Y with V, we obtain 



V v (x, y) = -|z^(x, y) - ^JT(y) . (48) 

Proceeding as before, with the help of ( JT5|) . the stationary equation ( 116cl) in the 
continuum is 

cu 2 [Z(x + 5, y + e 2 Ri _i) + Z(x -£,y + e 2 R- 1:1 ) -Z(x-e,y + e 2 R- 1 ^ 1 )- 

Z(x + e,y + e 2 R hl )] + 7 [V(x + e, y + £ 2 i?i,i) + V(x - e, y + e 2 i?_ li _i)+ (49) 

V(x + 5, y + £ 2 iVi) + V(x - s, y + e 2 ^-!,!) - 4V(x, y)] = . 

The leading order contribution is of 0(e 3 ), 

2cu 2 

V a;:c (x,y) = Z xy (x,y) . (50) 

7 

By integrating in x, we obtain 

2a; 2 

V x (x,j/) = Z y {x,y) + g(y) , (51) 

7 

where {?(y) is a suitable integration constant that will be determined by imposing the 
boundary conditions. Now, taking the derivative of (1181) w.r.t. x, and the derivative 
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of (1511) w.r.t. y, and summing the results, we obtain a differential equation for the 
behaviour of Z in the bulk, 

Z>xxxx(x, U) - ^rZyyix, V) = -Qy{y) ■ (52) 

This is the general equation, whose solution yields the behavior of the various fields in 
the bulk. 



4-4- Boundary conditions 

In this section we impose all boundary conditions. The various constraints allow not 
only to uniquely determine the behaviour in the bulk, but also to establish a link with 
the physically relevant observables, such as the temperature profile. Analogously to the 
previous section, we proceed into two steps by separately analysing the implications of 
( gfifl and of fTT6cj) . 

By setting i = j in ( |39i) . we obtain 

Yi f i — Y^j-i + Vi+i^ — Vj 5 j + 7 [— 2Zj + i i j + 2Zj i j + Zj + i ) j_i — Zj 5 j_i+ 

— Zj^+i] = . 

We recall that in order to avoid the complication of the absolute value in the denominator 
of (1281) . we have assumed that % > j. Accordingly, the use of (128!) requires to consider 
always the lower (by convention) triangle of all the matrices. In order to satisfy this 
condition, we exploit the antisymmetry of Z to obtain 

Y^j — Yj^-i + Vj +lji — Vj ; j + 7 (— Zj + x,i + — Zj^-i) = , 

and, in field variables, 

2<%) - Y(e, y + e 2 i?i,_i) + V(e, y + e^i.i) - T(y)+ 

7 (-Z(e, ?/ + e 2 R 1A ) + Z(2e, j/ + e 2 R 2fi ) - Z(e, y + e 2 i?i,-i)) = . 1 j 

The leading contribution of ( |54l) is 0(1), as expected for the diagonal terms, yielding a 
boundary condition for f2: 

n(y) = 2<S>(y)-T(y)=0 . (55) 

This last expression is just a local version of the virial theorem for the harmonic 
oscillators. 

The reader can verify that the leading term of PU|) in the upper diagonal (i = j — 1) 
does not give further information. However, adding the equation for the upper diagonal 
to (|53|) . we obtain, in field variables, 

Y(x + e,y + £ 2 R ltl ) -Y(x + e,y + £ 2 iV0 + T (v + ^0,2) - T(y)+ 

7 [-2Z(x + e, y + e 2 R hl ) -Z(x + e,y + e 2 R 1 - 3 ) + Z(x + 2e,y + e 2 R 2 , 2 )+ (56) 

Z(x + 2e,y + e 2 R 2fi ) -Z{x + e,y + e 2 R 1 ^ 1 )} = . 

This equation gives rise to two relations of leading order. The first is redundant as it 
confirms that Z is zero along the diagonal. The second relation is, instead, a differential 
equation for T(y), 

T y (y) + 1 Z xx (0 J y) = . (57) 
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It allows determining the temperature profile, once Z(x,y) has been determined. 
We now turn out attention to (116 c\i . Along the diagonal (i = j), it is 2 

T (2V M - Vi_i,i_i - Vi+i.i+0 + 2^ 2 (Z iil _ 1 - Zi +1|i ) = . (58) 

It is straightforward to show that the above equation is equivalent to 

|(V M - + u; 2 Z M _a = - J (59) 

where the integration constant J is nothing but the average heat flux. In fact, the 
energy flux Jj between the particles % — 1 and i is the sum of two contributions, a 
deterministic one jf^ , due to the interaction with the neighbours, and a stochastic one 
J- s \ originating from the collisions, 

Ji = Jf + J[ s) (60) 

with 

j\ d) = LU 2 (q^m) = u 2 Zi-i,i , (61) 

J1 S) = 1({PU)-(PD) = \ (Vi-1,-1 - V M ) , (62) 

where in both definitions we have adopted the convention that a positive flux corresponds 
to energy travelling from smaller to larger % coordinates. Accordingly, (|59|) states the 
well known physical fact that the heat flux is constant along the chain (i.e., independent 
of %). In the continuum limit, equation fl59|) writes 

| [T(y) - T(y + e 2 R ,. 2 )] + oo 2 Z(e, y + ^X-i) = ~J , (63) 

The leading contribution of the l.h.s. is of 0(e) and so must be J ( J = Je). As a result, 
we can write, 

u 2 Z x (0,y) = -J , (64) 

This is a relevant piece of information that will allow us to uniquely determine Z(x, y) 
in the bulk. 

Finally, for the upper diagonal (i = j + 1), (I16cl) becomes 

oo 2 [Z(0, y) - Z(2e, y + e 2 R 2 , 2 ) + Z(2e, y + e 2 R 2fl ) - Z(0, y + e 2 R , 2 )\ + 

7 [V(2e, y + e 2 ^) + V(2e, y + e 2 i? 2 , 2 ) - 2V(e, j/ + e 2 R 1A )] = , 1 j 

from where we obtain to leading order 

2^ 2 Z,(0,y)+ 7 V x .(0,y)=0, (66) 
that, by virtue of ( 1331) . implies 

V,(0,y) = 0. (67) 

For f ll6cl) . combinations of the diagonal and upper diagonal relations give no further 
information. 

| The reader can easily check that if one identifies T(— 1) with the left temperature T + and T(+l) 
with the right temperature T~, then the boundary terms in (|16cj) cancel to each other, namely 
(RV + VR) = 2k B T (R + 77S). 
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4-5. Solution of the equations 

In this section we solve the differential equations of covariance matrices, to leading order 
in e. From this solution we derive analytical expressions for the temperature profile and 
the energy flux. We start noticing that the function Q{y) appearing in (15TT) is identically 
equal to zero. This is seen by setting x = and using (1331) and (IBTj) . As a result, (l52l) 
simplifies to, 

4c; 2 

%xxxx(x, y) - —Z yy (x, y) = . (68) 

The form of the above equation suggests to look for a solution by the method of 
separation of variables. Furthermore, the numerical solution of the stationary solution 
fU6ap suggests that Z(x, y) = at the boundaries of the domain V. Therefore, we 
assume the following Ansatz 

Z(x, y) = Bn(x) sin [(3 n (y + 1)] , Pn=™. (69) 



B n . (70) 



which, upon substitution into ( 1681) . gives 
d 4 B n ( nnuj \ 2 

dx 4 

This is readily solved by finding the four roots of the associated characteristic 
polynomial. Two of the four eigenvalues having a positive real part would lead to 
an unphysical divergence in x and have to be discarded. Another constraint is imposed, 
by recalling that Z(0,y) = 0. Altogether, the coefficients B n can be written as, 

( ' nix<jj \ 

B n (x) = A n exp(-a n x) sin(a n x) , a n = ( ^— J . (71) 
Finally, the constants A n can be determined by imposing (164"]) 

Z(x, y) = } — — exp(-a n x) sin(a n a;) sin(/3 n (y + 1)) . (72) 

odd n 

The only remaining unknown, J7", can be finally determined by imposing that the 
temperature profile interpolates between T + and T_. By integrating (137)1 in y, we find 

T(y)=T- 7 [ Z xx (0,s) ds , (73) 



where we have identified T(0) = T = (T + + T )/2. By substituting expression (I72|) into 
( [731) and performing the integral term by term, we obtain 

T(y) = T + M. % oos(^(y + 1)) • (74) 

odd n Pn 

The value of J is obtained by imposing T(— 1) = T + . From ( [741) . it follows that 

v 7 odd n 
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£ "" 3/2 = ^ (|) , (76) 



V8 \2, 

odd n 

where the Riemann ^-function has been introduced, we obtain for J 

7T 3 CU 3 \ 1/2 AT 



J V 7 J 8( v / 8-l)C(3/2) ' (77) 
which corresponds to expression f[2"Uj) for the heat flux in the thermodynamic limit. 
Moreover, by substituting J into (174|) . we obtain expression (|T9|) for the temperature 
profile. Finally, the equation (ITT)) allows a unique identification of Z. From (1T2|) we find 
AT 

Z(x, y) = V n~ 3/2 e- anX sm(a n x) sin( f3Jy + 1)) , (78) 



5. Discussion and open problems 

Several comments are in order about the analytical results derived in the previous 
section, starting from the expression (120!) for the leading term of the heat flux. First 
of all, we see that the flux J is proportional to the temperature difference AT. This 
feature does not only apply to the leading term, but is a general property which follows 
from the harmonic nature of the underlying dynamics. In more general contexts, we 
expect a nonlinear response regime to exist. 

Moreover, J is independent of the strength of the coupling with the baths A. This 
can be physically understood by realizing that A plays the role of the inverse of a contact 
resistance. In the thermodynamic limit, the overall thermal resistance is the sum of the 
contact plus the bulk contribution which eventually dominates, no matter how small is 
A. Only for A = 0, the asymptotic regime cannot be attained (the system is isolated). 
The coupling A will presumably manifest itself when accounting for higher order terms. 

It is interesting to notice the inverse square root dependence of J on the rate 7 of 
internal collisions. The limiting values 7 — > and 7 — > 00, signal a crossover towards 
a regime characterized by a slower (faster) decay of J, respectively. This is the case, 
because 7 = corresponds to an integrable dynamics, while for 7 = 00 the decay of the 
heat flux is determined by higher order terms. 

As for the temperature profile, we should stress that the equation ( |T9l) represents 
the first example of an analytic expression obtained in the presence of anomalous heat 
conduction. This is all the way more important, by recalling that, as noticed in [20J, 
the temperature profile of this stochastic model is quite similar to that found in a 
purely deterministic system such as the purely quartic Fermi-Pasta-Ulam chain. Even 
more remarkably, T(y) is a parameter-free function. Indeed, once the profile is shifted 
around the average temperature and the temperature difference is rescaled to unity, 
the resulting shape Q(y) is not only independent of A but also of 7 and u>. This 
suggests that the temperature profile might be universal (at least in the limit of small 
temperature differences in truly nonlinear systems). Unfortunately pure numerics alone 




Figure 2. Temperature profile T(y) as given by the analytical expression (fT5| . for 
T + = 1.1, T_ = 0.9, lo = X = 7 = 1 (dashed curve). The solid curves correspond 
to the profile T nurn , obtained from the numerical solution of equations (|16aH16c"|) for 
N = 100, 200, 400, 800 . The finite-size deviations from (JT5J), 6T = T— T num , rescaled 
by TV 1 / 3 , are shown in the inset. 



is not sufficient to clarify this issue. Finally, we wish to comment on the singularity 
observed at the two extrema, namely for y close to —1 and 1. From ( fi~9i) . we find that 
for y = — 1 + Sy 

l/Sy 

5Q(y) « ~ ^ ' ( 79 ) 

odd n 

where the cosine has been approximated with a parabola and the sum has been limited 
to n < l/Sy, to prevent that the argument of the cosine becomes larger than 0(1). 
Altogether the above equation tells us that the profile is characterized by a square root 
singularity. 

Although our analysis has allowed us to determine an exact expression for the field 
2i(x,y) at leading order in the bulk, and thus for the temperature profile and the heat 
current in the steady state, the determination of the other fields V(x,y) and Y(x,y), 
requires the knowledge of higher-order terms. Indeed, the integration constant F{y) in 
( |47|) that helps determining V(x, y) and Y(x, y) cannot be obtained from our analysis. 
A comparison with numerics [23] suggests that Tiy) = 0, but none of the equations we 
have analysed in the previous section supports this observation. Presumably one should 
consider some other combinations of the equations f ll6aH16cl) . but the investigation is 
hindered by the fact that we are not entitled to use any information on the behaviour 
of higher-order terms. 

As a matter of fact, the estimation of the higher-order terms, starting from the 
leading corrections is a highly nontrivial problem, since such terms are likely to be 
non-analytic in the smallness parameter e. This is seen by comparing the analytical 
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Figure 3. Finite size deviation of the heat flux SJ = J — J n um, rescaled by TV 1 / 2 , 
as a function of the size of the chain N, for 7=1 (circles), 2 (triangles), 5 (pluses) 
and 10 (stars), other parameters as in the previous figure. The lines correspond to 
power-law fits, from which we extract that the corrective terms scale as — 5 J ~ N~ a 
with a = 0.927, 0.914, 0.939 and 0.947, respectively. 

results and the numerical solutions for finite chains. The first evidence is presented in 
figure El where we have plotted the analytical profile T(y) and the numerical ones T num 
computed for chains of different lengths N. From the data in the inset, we deduce that 
T — T num is approximately proportional to A" 1 / 3 . While this confirms the correctness of 
expression (fl~9]) . it also indicates that the leading correction is of 0(e 2 / 3 ). Nonanalytic 
corrections affect also the heat current. This is illustrated in figure [31 where we plot 
the difference between the numerical values J nU m and the leading-order term J, formula 
( f20|) . for different system sizes A and for various 7 values. In all cases, we see a clean 
power-law convergence to zero but the value of the power is systematically smaller than 
1, meaning once again that non- analytic higher-order terms in e exist. An appropriate 
scheme for the treatment of higher-order terms remains an open question. 
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